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Abstract. There is a class of statistical problems that arises in several 
contexts, the Lattice QCD problem of particle physics being one that 
has attracted the most attention. In essence, the problem boils down to 
the estimation of an infinite number of parameters from a finite number 
of equations, each equation being an infinite sum of exponential func- 
tions. By introducing a latent parameter into the QCD system, we are 
able to identify a pattern which tantamounts to reducing the system to 
a telescopic series. A statistical model is then endowed on the series, and 
inference about the unknown parameters done via a Bayesian approach. 
A computationally intensive Markov Chain Monte Carlo (MCMC) algo- 
rithm is invoked to implement the approach. The algorithm shares some 
parallels with that used in the particle Kalman filter. The approach is 
validated against simulated as well as data generated by a physics code 
pertaining to the quark masses of protons. The value of our approach 
is that we are now able to answer questions that could not be readily 
answered using some standard approaches in particle physics. 

The structure of the Lattice QCD equations is not unique to physics. 
Such architectures also appear in mathematical biology, nuclear magnetic 
imaging, network analysis, ultracentrifuge, and a host of other relaxation 
and time decay phenomena. Thus, the methodology of this paper should 
have an appeal that transcends the Lattice QCD scenario which moti- 
vated us. 

The purpose of this paper is twofold. One is to draw attention to a class 
of problems in statistical estimation that has a broad appeal in science and 
engineering. The second is to outline some essentials of particle physics 
that give birth to the kind of problems considered here. It is because of the 
latter that the first few sections of this paper are devoted to an overview 
of particle physics, with the hope that more statisticians will be inspired 
to work in one of the most fundamental areas of scientific inquiry. 
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1. INTRODUCTION AND OVERVIEW 

Lattice Quantum Chromodynamics, or Lattice 
QCD, is an actively researched topic in particle phy- 
sics. Many investigators in this field have received 
the Physics Nobel Prize, the 2004 prize going to 
Gross, Politzer and Wilczek, developers of the no- 
tion of "asymptotic freedom" that characterizes 
QCD. Underlying the Lattice QCD equations are 
issues of parameter estimation that have proved to 
be challenging. Essentially, one needs to estimate an 
infinite number of parameters from a finite number 
of equations, each equation being an infinite sum of 
exponential functions. 

The approach proposed here is Bayesian; it is dri- 
ven by a computationally intensive Markov Chain 
Monte Carlo (MCMC) implementation. However, to 
invoke this approach, we need to introduce a latent 
parameter and then explore the "anatomy" of the 
QCD equations. This reveals a pattern, which when 
harnessed with some reasonable statistical assump- 
tions provided a pathway to a solution. The inferen- 
ces provided by our approach were successfully vali- 
dated against simulated as well as real data. However, 
the real value of our approach is that it is able to ans- 
wer questions that could not be answered using some 
of the conventional approaches of particle physics. 
The approach can therefore be seen as an addition 
to the lattice field theorists' data analysis tool kit. 

The structure of the Lattice QCD equations is 
not as specialized as one is inclined to suppose. In- 
deed, such equations also appear in other contexts of 
engineering, physics, nuclear magnetic imaging and 
mathematical biology where they go under the la- 
bel of "exponential peeling;" see Section 3.1. Our 
focus on the physics scenario is due to the fact that 
this is how we got exposed to the general problem 
addressed here. 

This paper is directed toward both statisticians 
and physicists, and could serve as an example of 
the interplay between the two disciplines. The for- 
mer may gain an added appreciation of problems in 
modern physics that can be addressed via statistical 
methods. In the sequel, they may also get to know 
more about particle physics and the beautiful theo- 
ries about it that Mother Nature has revealed. It is, 
with the above in mind, that Section 2 is devoted to 
an overview of aspects of particle physics, its asso- 
ciated terminology and the awe inspiring discover- 
ies about it. Reciprocally, the physicists may benefit 
by exposure to some modern statistical technologies 
that can be brought to bear for addressing problems 
that may have caused them some consternation. 



Section 2 gives an overview of some essentials of 
particle physics, and the ensuing Lattice QCD equa- 
tions. This section, written by a nonphysicist (NDS) 
but reviewed by a physicist (FXL), has been devel- 
oped by fusing material from a variety of sources, 
some notable ones being Pagels (1982), Dzierba, 
Meyer and Swanson (2000), Yam (1993), Riordan 
and Zajc (2006) and Frank Wilczek's (2005) Nobel 
lecture. Interjected throughout this section are a few 
comments of historical interest; their purpose is to 
inform a nonphysicist reader about the individuals 
who have contributed to the building of a magnif- 
icent edifice. Section 2 concludes with a graphical 
display of the structure of matter via a template 
that is familiar to statisticians, in particular, those 
working in network theory and in reliability. 

Section 3 pertains to an anatomy of the Lattice 
QCD equations and the resulting mathematical pat- 
tern that it spawns. It is not necessary to read Sec- 
tion 2 (save perhaps for an inspection of Figure 5) 
in order to read Section 3, which is where this paper 
really begins; indeed, Section 2 could have been del- 
egated to an Appendix. Section 3 is a foundation for 
the rest of the paper. It is here that the inferential 
problem is introduced along with its accompanying 
notation and terminology. Section 3.1 gives a broad 
overview of the several other scenarios in science and 
engineering where the Lattice QCD type equations 
also arise. Of particular note are the several exam- 
ples in mathematical biology wherein the QCD like 
equations are often discussed. 

Section 4 pertains to the statistical model that the 
material of Section 3 creates, and an outline of the 
MCMC approach that is used to estimate the para- 
meters of the model. These are the parameters that 
are of interest to physicists and other scientists. Sec- 
tion 5 pertains to validation against simulated and 
actual data and proof of principles. Section 6 per- 
tains to some suggestions for extending the work do- 
ne here, and strategies for overcoming some of the en- 
countered difficulties. Section 7 concludes the paper. 

Since the Lattice QCD equations can be seen as 
a prototype for similar equations that arise in other 
scientific endeavors, this paper also serves as an in- 
vitation to other statisticians to develop approaches 
for solving such equations using methods more so- 
phisticated and/or alternate to the one we have en- 
tertained. 

2. ESSENTIALS OF PARTICLE PHYSICS 

The smallest quantity of anything we can see or 
feel is a molecule, and all matter is made up of mo- 
lecules, which in turn are made up of atoms. Molecu- 
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Fig. 1. Architecture of a carbon atom. 

les and atoms are called particles, and the physics 
that describes the interactions between the parti- 
cles is known as particle physics; see, for example, 
Griffiths (1987). 

An atom consists of electrons, which carry a neg- 
ative charge, and the electrons are centered around 
a nucleus that is made up of protons that carry a po- 
sitive charge, and neutrons that carry no charge. Fi- 
gure 1 illustrates the architecture of a carbon atom 
which has six electrons, six protons and six neutrons; 
it is denoted g 2 C. 

The protons and the neutrons are held together 
within the nucleus by a nuclear glue called the pion. 
Similarly, the protons and the electrons are held to- 
gether within the atom by a glue called the pho- 



ton. The pions are said to be carriers (or media- 
tors) of the strong force (or the nuclear force), and 
the photons are carriers of the electromagnetic force. 
Physicists look at the nuclear glues as force carrying 
particles, and thus collectively regard the electrons, 
the neutrons, the photons, the pions and the pro- 
tons as subatomic particles. Figure 2 displays the 
structure of matter as understood around the 1946 
time frame. The dotted lines of Figure 2 indicate the 
glued members. 

In 1911, when Rutherford announced the struc- 
ture of the atom, the existence of electrons and pro- 
tons was known. The neutron, as a major constituent 
of the nucleus, was discovered in 1932 by Chadwick, 
and the pion was discovered in 1946. But these dis- 
coveries were just the tip of the iceberg. Many more 
subatomic particles have subsequently been discov- 
ered. Collectively, these subatomic particles are now 
called hadrons. Physicists speculate that there exist 
an infinite number of such hadrons. This discovery 
of hadrons was made possible by accelerators, which 
are essentially microscopes for matter. 

The invention of the accelerators opened up the 
subnuclear world with the experimental discovery of 
thousands of new particles. The question thus arose 
as to what the hadrons could be saying about the 
ultimate structure of matter. 

2.1 The Quark Structure of Matter 

The current view is that hadrons are composite 
objects made out of more fundamental particles cal- 
led quarks, and no one has ever seen a quark! This 
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Fig. 2. The structure of matter (circa 1946). 
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Orbit of Quark 2 



Fig. 3. Illustration of a quark orbit. 



point of view came about in the early 1960s when 
Murray Gell-Mann discovered that the hadrons or- 
ganized themselves into classes (or families) based 
on a mathematical symmetry. An easy way to under- 
stand why this organizational principle worked is to 
assume that the hadrons are made up of quarks, only 
three of which were needed to build the hadrons. 
These quarks were named the up quark, the down 
quark and the strange quark. For example, a proton 
has two up quarks and one down quark, whereas 
a neutron has two down quarks and one up quark. 
In general, every hadron is made up of quarks that 
orbit around each other in a specific configuration, 
each configuration resulting in a hadron. Figure 3 is 
an illustration of a quark orbit. 

Since there could be several orbit configurations, 
there ought to be an infinite number of hadrons. 



The essence of Gell-Mann's idea is that hadrons are 
bound states of quarks, just like how the atoms 
are bound states of electrons, neutrons and pro- 
tons. Furthermore, Gell-Mann postulated that there 
ought to exist a force carrying particle, called the 
gluon, that holds the quarks together. The gluon is 
said to be the carrier of the strong force. Figure 4 
illustrates the quark structure of a hadron. 

The quark model was purely a theoretical cons- 
truct. Its validity was affirmed when Gell-Mann used 
it to postulate in 1962 the existence of a particle 
never seen before. This was a scientific breakthrough 
of the highest order! It showed that discoveries in 
physics can come from mathematical patterns — not 
just the laboratory. For unraveling the mathematical 
symmetries of the hadron, Gell-Mann received the 
1969 Nobel Prize in Physics. 
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Fig. 5. Matter as a coherent system. 



Figure 5 gives a pictorial representation of the 
quark structure of matter using a template that is 
familiar to statisticians. It represents an atom as 
a coherent (or logical) system with quarks as the 
basic building blocks of the system. The logic sym- 
bols of "and" and "or" are represented by and re- 
spectively. The neutrons and the protons can be re- 
garded as subsystems, and the gluons, photons and 
the pions that link the quarks, the nucleus and the 
electrons can be seen as the structure (or link) func- 
tions of the system (cf. Barlow and Proschan, 1975). 
These are the carriers of the strong force and the 
electromagnetic force, respectively. Figure 5 contains 
Gell-Mann's famous quote that "everything that is 
not forbidden is compulsory;" the logical systems 
analogue to this quote is the notion of "irrelevance." 

2.2 Quantum Chromodynamics and Lattice QCD 

The theory of QCD can be thought of as a recipe 
for producing hadrons from quarks and gluons. Since 
quarks and gluons make up most of the known mass 
of the physical world, unraveling the quark structure 
of matter is the key to an understanding of the phys- 
ical world, and thus the importance of the subject 
of this paper. 

The QCD theory was successful in enunciating the 
properties of the hadrons. However, its complexity 
made its use for predicting unobservable quantum 
quantities, like quark masses, almost impossible. This 



is because solving the QCD equation (which is just 
one line) by analytical methods is difficult. The cur- 
rent approach is to solve the QCD equation numer- 
ically, by discretizing it over a space-time lattice. 
Lattice QCD refers to the representation of space- 
time scaffold in four dimensions wherein the 
quarks rest on the connecting sites, and the gluons 
as connections between the lattice points. 

The scaffold is first restricted to a finite volume; 
it is then replicated with periodic boundary condi- 
tions. All this entails on the order of 100 million 
billion arithmetic operations on typical lattices; this 
is one example as to why physicists need supercom- 
puters. Lattice QCD has been able to explain as to 
why a free quark has not been seen and will not be 
seen; this is because it will take an infinite amount 
of energy to isolate a quark. 

Lattice QCD, being an approximation to the QCD, 
improves as the lattice points increase indefinitely 
and as the volume of the lattice grid expands. In 
so doing it opens up avenues for statistical methods 
to enter the picture. Physicists have explored some 
of these avenues, one of which is the focus of this 
paper; see Section 3 below. 



THE UNDERLYING PROBLEM: 
EQUATIONS 



QCD 



With Lattice QCD, an archetypal scenario is the 
estimation of an infinite number of parameters from 
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a finite number of equations. The ieft-hand side of 
each equation is the result of a physics based Monte 
Carlo run, each run taking a long time to complete. 
Thus, there are only a finite number of runs. For 
example, a meson correlator, G(t\-), takes the form 
(cf. Lepage et al., 2002) 

oo 

(3.1) G(t\-)=J2A n e- E » t for t = 0,1, 2,..., 

n=l 

where the parameters A n denote the amplitude, and 
E n denote the energy. Also, E\ < E2 < ■ ■ • <E n <- ■ ■. 

Interest centers around the estimation of A n 
and E n , n = 1,2, ... , based on G(t\-), estimated as 
G(t\-),t = 0, 1, . . . , k, for some finite k [23 in the case 
of Lepage et al. (2002)]. The physics codes which 
generate the G(t|-)'s do not involve the A n 7 s and 
the En 8, and are autocorrelated, thus the label "cor- 
relator." The physics codes also provide estimates of 
the autocorrelation matrix. 

Deterministic approaches to solve for the A n 's and 
the E^s cannot be invoked, and statistical approa- 
ches involving curve fitting by chi-square, maximum 
likelihood and empirical Bayes have proved to be un- 
satisfactory (cf. Morningstar, 2002). For an apprecia- 
tion of these efforts, see Lepage et al. (2002), Fiebig 
(2002) and Chen et al. (2004); the latter authors 
propose what they call a "sequential empirical Bayes 
approach." However, empirical Bayes approaches use 
observed data to influence the choice of priors, and, 
as asserted by Morningstar (2002), are a violation 
of the Bayesian philosophy. Indeed, Fiebig (2002) 
states that "Bayesian inference has too long been 
ignored by the lattice community as an analysis 
tool. . . . The method should be given serious con- 
sideration as an alternative for conventional ways." 

Bayesian approaches alternate to ours have been 
considered by Nakahara, Asakawa and Hatsuda 
(1999). These authors entertain the use of maxi- 
mum entropy priors, but, as claimed by Lepage et al. 
(2002), the accuracy of their estimator of Ei is in- 
ferior to those obtained using other approaches. Be- 
cause priors based on the principle of maximum en- 
tropy result in default priors, such priors also violate 
the Bayesian philosophy. The approach of Lepage 
et al. (2002) is Bayesian in the sense that prior in- 
formation is used to augment a chi-square statistic 
which is then minimized. We find this work valuable 
because it articulates the underlying issues and pro- 
vides a framework for examining the anatomy of the 
QCD equations, which enables us to identify a pat- 
tern, which in turn enables us to invoke the Bayesian 
approach we propose. 



3.1 Relevance to Other Scenarios in Science and 
Engineering 

The Lattice QCD architecture of equation (3.1) is 
not unique to physics. They occur in several other 
scenarios in the physical, the chemical, the engi- 
neering and the biological sciences, a few of which 
are highlighted below. Most attempts at estimation 
of the underlying parameters have involved least 
squares or numerical techniques based on local lin- 
earization with iterative improvements. Besides lack- 
ing a theoretical foundation vis-a-vis the require- 
ment of coherence (cf. Bernardo and Smith, 1994, 
page 23), techniques have proved notoriously unre- 
liable and not robust to slight changes in the exper- 
imental data (cf. Hildebrand, 1956). 

Mathematical biology: exponential peeling in com- 
partment systems When considering radioactive tra- 
cers used for studying transfer rate of substances in 
living systems (cf. Robertson, 1957; Rubinow, 1975, 
page 125), sums of exponentials are encountered. 
Here, the G(t\-) of equation (3.1) represents the con- 
centration of a substance, the t's are integer values 
of time, and the A^s and the E^s are constants that 
need to be estimated. Here interest generally centers 
around the case of n — 2, and the coefficients A n 
and E n of equation (3.1) are negative. An ad hoc 
graphical procedure called the method of exponential 
peeling is used to estimate the parameters (cf. Smith 
and Morales, 1944, Perl, 1960; van Liew, 1967). 

Some other scenarios in biology where the Lattice 
QCD type equations appear are in bone metabolism 
studies and cerebral blood flow (cf. Glass and de 
Garreta, 1967), and in biological decay (cf. Foss, 
1969). In the latter context, Dyson and Isenberg 
(1971) consider for fluorescence decay an equation 
of the type 

m 

j'=l 

where y{t) represents "moments of the fluorescence," 
ctj's the amplitudes [the A n 's of equation (3.1)], 
and the Tj's are time constants corresponding to 
the E n 's of equation (3.1). Here the ay's are zero 
for j > m + 1 . 

Gene expression data When considering a time 
series of gene expression data (cf. Giurcaneanu et al., 
2005), a system of equations paralleling that of equa- 
tion (3.1) arises again. In this context G(t\-) repre- 
sents "mRNA concentrations" as a function of time, 
and the parameters A n and E n describe interactions 
between the genes. In the gene expression context, 
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as in the Lattice QCD context, the parameters E n 
are increasing in n. 

Nuclear magnetic resonance (NMR) NMR exper- 
iments often generate data that are modeled as the 
sum of exponentials (cf. Bretthorst et al., 2005). 
Experiments relying on NMR to probe reaction ki- 
neticis, diffusion, molecular dynamics and xenobi- 
otic metabolism are some of the applications where 
parameter estimates provide insight into chemical 
and biological processes. See, for example, Paluszny 
et al. (2008/09) who study brain tissue segmentation 
from NMR data. 

Here one considers equations of the type 
m 

di = C + ^2 Aj exp{-a>jti} + rij, 
i=i 

where m is the number of exponentials and di a data 
value sampled at tj. The parameters of interest are 
the decay rate constants ay, the amplitudes Aj and 
the constant offset C. The n s 's are the error terms. 

Electromechanical oscillations in power systems 
Equations entailing the sum of exponentials are also 
encountered in the context of low frequency electro- 
mechanical oscillations of interconnected power sys- 
tems, the impulse response of linear systems in net- 
works, ultracentrifuge and a host of other relaxation 
and time-decay phenomena (cf. Dyson and Isenberg, 
1971). For example, in the electromagnetic oscilla- 
tions scenario, Sanchez-Gasca and Chow (1999) en- 
counter an equation analogous to our equation (3.1) 
with G(t\-) denoting a signal and A n connoting a sig- 
nal residue associated with the "mode" E n . 

To summarize, the relationships of the type given 
by equation (3.1) arise in so many contexts of science 
and engineering that it seems to be quintessential, 
and almost some kind of law of nature. The Lat- 
tice QCD problem considered here can therefore be 
seen as a prototype and a convenient platform to ex- 
posit a statistical problem of general applicability. In 
most of the application scenarios described above, 
statistical methods have been used, many ad hoc, 
some empirical Bayesian and a few Bayesian (un- 
der the rubric of maximum entropy) . Many of these 
methods have not exploited an underlying telescopic 
pattern in these equations which makes an appear- 
ance when a latent parameter is introduced into the 
system, and inference about the latent parameter 
made. 

3.2 Anatomy of the Lattice QCD equations 

An examination of equation (3.1) yields the fol- 
lowing boundary conditions. G(0|-) = Yl^Li A n , im- 



Values of t f 




A, A 2 A 3 



Number of^4 w 's 
Fig. 6. Number of A n 's as a function oft. 

plying that the A^s are constrained. When t — > oo, 
G(t\-) = 0, which implies that for large values of t, A n 
and the E n cannot be individually estimated. Thus, 
simulating G(t\-) for large t does not have a payback; 
consequently, it is futile to do such a simulation. 

Since the E^s increase with n, we may, as a start, 
reparameterize the E^s as E n — E n —\ = c, for some 

unknown c, c> 0, for n = 2, 3, It will be argued 

later, in Section 6.1, that c is a latent parameter. 
Thus, 

(3.2) E n = E x + (n - l)c, n = 2,3,... 

with E\ and c unknown. With the above assumption 
in place, a parsimonious version of the Lattice QCD 
equation takes the form 

oo 

(3.3) G(t\-)=e- Elt J2 A ne~ {n ~ 1)ct , 

n=l 

£ = 0,1,2,.... 

With c fixed, the parsimonious model given above 
reveals the following features: 

(a) When t is small, the number of A^s entering 
equation (3.3) is large; indeed, infinite when t = 0. 

(b) When t is large, the number of ^4 n 's we need to 
consider is small, because the combination of a large t 
with any n will make the term A n exp(—(n — l)ct) 
get small enough to be ignored. 

(c) Moderate values of t and n will also make the 
above term small, causing A n to be irrelevant. 

Figure 6 illustrates the feature that as t gets large, 
the number of A^s one needs to consider gets small. 

As a consequence of the above, for any fixed c, we 
can find a t\ such that in the expression 

e -tEi[ Al+A2e -ct + A3e -2ct + ... 

+ A n e-^ ct + •••], 
all the terms, save for A±, are essentially zero. 
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Similarly, we can find a t 2 , t 2 < t±, such that all 
the terms save for A\ and A2e~ ct2 get annihilated. 
Continuing in this vein, there exists a sequence t k < 
tk-i < ■ ■ ■ < t% < t±, such that all that matters are 
the terms associated with A± , A 2 , ■ ■ ■ , A^ • In what 
follows, we suppose that k is specified. 

Thus, for any fixed c and k, with t\ > t 2 > • • • > tk 
chosen in the manner described above, our parsimo- 
nious version of the Lattice QCD equations telescope 
as follows: 



(3.4) 



G(t 2 \-) 
G(t k \-) 



-EltlA 



e- E ^ 2 (A 1 + A 2 e- ct2 ), 



-E\t n 



(A 1 + A 2 e~ ctk + 



+ A k e- {k ~ l)ctk ). 

To summarize, by introducing the constant c, fixing 
a k, and identifying an underlying pattern in the 
Lattice QCD equations, we have reduced the prob- 
lem to the case of k equations and (A; + 2) unknowns, 
A\, . . . , Ak,E% and c. The choice of what k to choose 
is determined by the number of physics code based 
estimates G(t),t = 0, 1, . . . , k, that can be done and 
are available. 

4. STATISTICAL MODEL: SOLVING THE QCD 
EQUATIONS 

Many have expressed the view that it would be 
considered good progress if trustworthy estimates of 
just Ai,A 2 ,E\ and E 2 can be had. The other pairs 
(^3, £3), (A4, £4), . . . , can be considered later; see 
Section 6. Thus, we start by focusing attention on 
the first two equalities of equation (3.4); that is, the 
case k = 2 and some fixed c. Specifically, we consider 

G(t 1 \A 1 ,E 1 ) = e~ Eltl A 1 and 

(4.1) 

G(t 2 \Ai ,E 1 ,A 2 ,c) = e~ Elt2 (A l + A 2 e 



If Hi = G(ti\-),i = 1, 2, denotes the physics code based 
evaluations of G(U \ ■) , then our aim is to estimate A\ , 
E±,A 2 and c, in light of y\ and y 2 . To set up our like- 
lihoods, we take a lead from what has been done 
by Nakahara, Asakawa and Hatsuda (1999), and by 
Lepage et al. (2002), to write 

Y x = G(tx\-) +£1 and 

Y 2 = G(t 2 \-)+e 2 , 
where £j ~ N(0, crf),i = 1, 2, and Corr(ei, e 2 ) = p\ 2 . 



Besides providing y\ and y 2 , the physics codes also 
provide a\ , a\ and p\ 2 . As a consequence, the sta- 
tistical model boils down to the bivariate normal 
distribution, 



(4.3) 



~Y{ 




'G(h)' 




T2_ 
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Writing out a likelihood function for the unknowns 
A\,E\,A 2 and c, based on equation (4.3), is a straight- 
forward matter. However, we need to bear in mind 
that since the parameters A\ and E\ appear in both 
G(ti\-) and G(t 2 \-), both y\ and y 2 provide informa- 
tion about A\ and E±, with y 2 providing informa- 
tion about A 2 and c as well. To exploit this feature, 
we construct our likelihoods based on the marginal 
distribution of Y\, and the conditional distribution 
of Y 2 given Y\. That is, on 

(4.4) Yi~ N(A ie ~ Eltl ,af) 

and 



°2 , 



(Y 2 \Y 1 =y 1 )~N(G(t 2 \-)+p 12 — (y 1 -G(t 1 \-)), 



^2 2 (1-P? 2 ) 



(4.5) 



Specifically, the likelihood of A\ and E\ , with y\ fi- 
xed, is 

fAR\ CIA TP \ \ fel-^"^) 2 ' 

(4.6) £(^i,£i;yi) oc exp ^ 

and the likelihood of A±, E±, A 2 and c, with y 2 fixed, 
and the effect of y\ incorporated via the posterior 
distribution of A\ and E\ , is of the form 

C(A 1 ,E 1 ,A 2 ,c;y 1 ,y 2 ) 



oc exp 



(4.7) 



{y 2 -(e- Elt2 (A 1+ A 2 e~ ct2 )) 



+ P12 — (yi - Aie 



[2cr 2 2 (l-p? 2 )]- 



In the above development, the covariance matrix is 
provided by the physics code. As suggested by a ref- 
eree, a deeper investigation of this matrix may be 
called for, because with increasing t, the variances 
are likely to increase, posing computational chal- 
lenges to the proposed approach. 

4.1 Specification of the Prior Distributions 

To implement our Bayesian approach, we need to 
make assumptions about conditional independence, 
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and assign prior distributions for the unknown pa- 
rameters. The priors that we end up choosing in Sec- 
tion 5 are not based on knowledge of the underlying 
physics, but are proper priors based on an appreci- 
ation of the material in Morningstar (2002), Lepage 
et al. (2002) and Fleming (2005). 

The Ai's are supposedly between and 1, and no 
relationship between them has been claimed. Thus, 
it is natural to assume that A\ and A<i are apriori 
independent, and have a beta distribution on (0, 1) 
with parameters (a, (3); we denote this as B(Af, a, f3), 
i = l,2. The relationship between E\ and c is less 
straightforward. We conjecture that the larger the E\ , 
the smaller the c, and that E\ can take values over 
(0, oo). It is therefore reasonable to assume that the 
prior on E\ is a gamma distribution with scale pa- 
rameter r\ and shape parameter A; we denote this 
by Q(Ei;rj,X). Some other meaningful choices for 
a prior on E\ could be a Weibull, or a Pareto, the lat- 
ter being noteworthy as a fat-tailed distribution. To 
encapsulate the dependence between E\ and c, we 
suppose that, given E\, c has a uniform distribution 
over (0,uj/Ei), for some oj > 0. Finally, we also as- 
sume that Ei and c are independent of all the A^s. 

The above choice of priors, with user specified hy- 
perparameters a, /3, oj, A and r/, is illustrative. In 
principle, any collection of meaningful priors can be 
used, since the ensuing inference is done numerically 
via a Markov chain Monte Carlo (MCMC) approach. 

Lepage et al. (2002), and also Morningstar (2002), 
seem to use independent Gaussian priors for the pa- 
rameters in question — see equations (8) and (11) 
respectively. Indeed, Morningstar (2002) makes the 
claim that "practitioners often restrict the choice of 
a prior to some familiar distributional form." The re- 
stricted parameter space makes the choice of Gaus- 
sian priors questionable. An overview of how the 
MCMC is invoked here is given next. 

4.2 An Outline of the MCMC Excercise 

The telescopic nature of the Lattice QCD equa- 
tions suggests that the MCMC will have to be con- 
ducted in the following three phases: 

Phase I. Using E^ as a starting value and y\ as 
data, obtain the posterior distribution of A\ and E\ 
via equation (4.6) as the likelihood, and 1,000 iter- 
ations of the MCMC run. 

Phase II. Using c*- -* as a starting value, and 1/2 &s 
data, obtain a sample from the posterior distribu- 
tion of Ai, E\, A2 and c via the likelihood of equa- 
tion (4.7), and 1,000 iterations of the MCMC run. 
Sample values of Ai and Ei from their posterior dis- 



tributions obtained in Phase I will serve as the priors 
of Ai and Ei in Phase II. Since the parameters Ai 
and Ei reappear in the likelihood of equation (4.7) 
as the mean of 7/2 , Phase II of the MCMC run cap- 
tures the effect of 1/2 on these parameters. The effect 
of yi was captured in Phase I. 

Phase III. Repeat Phase I and Phase II m times 
using new starting values of Ei and c to produce 
a sample of size m from the posterior distribution 
of Ai, Ei, A2 and c, with yi and y2 as the data. 

The MCMC exercise described above is routine, 
but computer intensive and entails 12 steps, six in 
each phase, and this too for a highly curtailed ver- 
sion of the Lattice QCD equations. The details of 
how this is done could be interesting, because they 
involve some discretization of the simulated poste- 
rior distributions, and working with individual sam- 
pled values reminiscent of that done in particle Kal- 
man filtering (cf. Gordon, Salmond and Smith, 1993). 
Thus, we label our approach as Particle MCMC. 
More details are given in Landon (2007), and the 
method illustrated in the Appendix. The software 
can be downloaded at http://www.gwu.edu/~stat/ 
irra/Lattice_QCD.htm. 

4.3 A Caveat of the Proposed Scheme 

The caveat mentioned here stems from the fea- 
tures that c has been fixed, and that the MCMC 
runs are centered around fixed values of yi and 2/2 • 
To see why, recall that our parsimonious version of 
the Lattice QCD equations [see equation (3.4)] is 
based on those ij's for which the exponential terms 
vanish; however, the ij's are determined by a fixed 
value of c. Thus, any change in the value of c will 
bring about a change in the values of ti, and, as 
a consequence, the Lattice QCD equations will also 
have to be different. This would be tantamount to 
obtaining new values of the y^s. However, all the 
likelihoods in the MCMC runs are based on fixed 
values of the yi's; see equations (4.6) and (4.7). But 
a change in the value of c is inevitable, because 
in Phase II of the MCMC run one iterates around 
sampled values from the posterior distribution of c, 
so that the initial systematically gets replaced 
by c« C ( 2 ),..., C (booo). 

A way to overcome this caveat is to recognize that 
for any > , i = 1, 2, . . . , the exponential terms 
mentioned above will continue to vanish, so that 
any specified values of yi will continue to satisfy the 
right-hand side of equation (4.2). 

A strategy to ensure that the successively gen- 
erated values of c^, i = 1,2,..., will tend to be 
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greater than is to pick small values of for 
each of the m iterations of Phase III of the MCMC 
algorithm. During the course of the MCMC runs, 
should one encounter a generated value of cW that 
is smaller than , then one should discard the 
so-generated value cW, and generate another value 
of cW . Hopefully, the number of discarded cW 's will 
not be excessive, but if they are, then the start- 
ing value c(°) should be decreased, and new values 
of t\ and ti obtained. This of course would be tan- 
tamount to obtaining new values of y\ and yi as 
well. 

5. PROOF OF PRINCIPLE: VALIDATION 
AGAINST DATA 

We first validate the accuracy of our approach 
against simulated data. For this, we choose A\ = 0.8, 
A 2 = 0.6, A 3 = 0.4, A A = 0.2, A 5 = 0.1, and Ai = 
for i > 6. We also choose E\ = 0.9 and c = 0.5. Using 
these values in equation (3.3), we compute G(t), for 
t = 1, 2, . . . , 12; these are shown in column 3 of Ta- 
ble 1. Since Y t = G(t\-) + e u with e t ~ 7V(0, of) [see 
equation (4.2)], we generate yi,...,yn, assuming 
the £f 's are independent, with at = 0.001 x G(t) x t; 
these are shown in column 4 of Table 1. We next 
identify those t's for which the leading exponen- 
tial terms vanish. These happen to be t\ at t = 12, 
t-2 at t = 6, and £3 at t = 4; see column 2 of Ta- 
ble 1. Our aim is to invoke the methods of Section 4 
on the entries of Table 1, to see if the constants 
specified above can be returned. With the above 
in place, Phases I, II and III of the MCMC run 
were made arbitrarily choosing the hyperparameters 
a = (3 = 77 = A = uj = 1, and m = 1,000. 



Table 1 

Simulated data for validating approach 



Time t 


Index ti 


G(t) 


Vi 


1 




0.54874373 


0.54900146 


2 




0.17764687 


0.17756522 


3 




0.06387622 


0.06373037 


4 


t 3 


0.02422158 


0.02414992 


5 




0.00945326 


0.00952723 


6 


t 2 


0.00375071 


0.00377058 


7 




0.00151265 


0.00151498 


8 




0.00060552 


0.00061147 


9 




0.00024486 


0.00024698 


10 




0.00009923 


0.00009821 


11 




0.00004026 


0.00004007 


12 


ti 


0.00001635 


0.00001625 



5.1 Results Based on Simulated Data 

Figure 7(a) and (b) shows the posterior distribu- 
tions of E\ based on j/12, and on y\2 and ye, re- 
spectively. Recall that j/12 corresponds to t\ , and ye 
corresponds to ti- Note that the posterior distribu- 
tion of Figure 7(a) becomes the prior distribution 
for the construction of the posterior distribution of 
Figure 7(b). Both the distributions of Figure 7 indi- 
cate a modal value of 0.9, suggesting a tendency to 
converge to the true value of E\. Furthermore, the 
difference between the two distributions is not very 
great, suggesting that ye may not be contributing 
much toward inference for E\ , beyond that provided 
by J/12- 

A similar feature is revealed by the posterior dis- 
tributions of A\, shown in Figures 8(a) and (b). 
These distributions have a modal value of 0.8, sug- 
gesting again a convergence to the true value of A± . 
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0.03 



0.025 



0.015 




0.005 



Fig. 8. Posterior distribution of Ai. (a.) Posterior of A\ based on j/12. (b) Posterior of Ai based on y\2 and y§. 



Figures 9 and 10 show the posterior distributions 
of A2 and c, based on yi2 and y§. Their modal 
values of 0.6 and 0.5 suggest convergence of the 
posteriors to their true values. Thus, based on this 
simulation exercise, we may claim that, despite an 
arbitrary choice of hyperparameters, the proposed 
MCMC procedure is able to show recovery of the 
input values of A±, E\, A2 and c to a meaningful 
degree of accuracy. 

5.1.1 Sensitivity of posteriors to priors In this sec- 
tion we explore the sensitivity of the posterior dis- 
tributions of Ai, Ei, A2 and c when the hyperpa- 
rameters of their prior distributions vary. We also 
explore the effect of using a thick-tailed prior distri- 
bution for Ei, in particular, a Pareto distribution, 
instead of the gamma distribution used before. 

Figure 11 shows the posterior distributions of Ei 
for different values of the scale A and shape rj pa- 
rameters of its gamma prior. Verify that the poste- 
rior distributions get centered around its true value 
of 0.9 even when the prior mean is as large as 10. 
The values of the chosen hyperparameters are indi- 
cated in the legend accompanying Figure 11. 

In Figures 12 and 13 we show the posterior dis- 
tributions of A\ and A2 for different values of the 
hyperparameters a and /3; see the legend accompa- 
nying these figures. Whereas the posterior distribu- 
tion of A2 appears to be very robust against the 
various choices for its prior distributions, the pos- 
terior distribution of Ai shows some sensitivity — 
albeit minor — to the choice of its priors. These pri- 
ors are centered at (in the case of A2) and around 
(in the case of A{) their true values of 0.6 and 0.8, 
respectively. 




0.1 0.2 0.3 0.4 0.5 06 07 08 09 1 

Fig. 9. Posterior of A2 based on yi2 and y$. 



0.12 




0.06 



0.01 



0.02 



01 02 03 04 05 06 0.7 0.8 0.9 1 

Fig. 10. Posterior of c based on yi2 and ya- 

Since the prior on c is a uniform on (0,uj/Ei), 
changing the value of ui would simply change the 
range of values that c can take. It will not change 
the shape of the posterior distribution of c. Finally, 
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Fig. 11. Posterior distribution of Ei with different values of rj and X. 



0.025 



0.02 



0.015 



0.01 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Fig. 12. Posterior distribution of Ai with different values of a and p. 



a use of the Pareto as a prior for E\ results in a pos- 
terior distribution that looks much like that of Fig- 
ure 7 produced by a gamma prior. This result — not 
illustrated here — is true irrespective of the choice 
of the hyperparameters of the Pareto prior. Indeed, 
the Pareto prior for E\ indicates a higher degree of 
robustness of its resulting posterior as compared to 
the gamma prior. 

Overall, it seems to be the case that the proposed 
procedure is robust to the choice of priors, and that 
the resulting posteriors converge to their correct val- 
ues no matter the choice of priors. 



5.2 Validation Against Physics Code Data 

In this section we validate our approach using data 
pertaining to a pion that has been generated by 
a physics based code. These data are given in Table 2 
and parallel those of Table 1, save for the fact that 
the data run from t = 2 to t = 13, and that G(t),t = 
2, . . . , 13, is not known. However, the Y^s and their 
associated errors are provided by the code, the errors 
being a proxy for a\. The choice of ti,*2 an d H is 
based on the following consideration. By default, t\ 
has to be the largest t for which the data are availab- 
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Table 2 
Physics code based data on a pion 



t 


Index ti 


Y t 


Errors 


2 




0.043865236 


0.00013635 


0.00014836 


3 




0.009347211 


0.00008205 


0.000089027 


4 


t 3 


0.00406969 


0.000051302 


0.000057832 


5 




0.002187666 


0.000031545 


0.000034867 


6 




0.001252858 


0.000018805 


0.000018559 


7 




0.000735911 


0.000011131 


0.00001124 


8 




0.0004358 


6.6393E-06 


6.8252E-06 


9 




0.00025829 


0.000004049 


4.3108E-06 


10 




0.000153161 


2.4808E-06 


2.6302E-06 


11 




9.1412E-05 


1.5264E-06 


0.000001683 


12 




5.552E-05 


9.5081E-07 


1.0741E-06 


13 


h 


3.54336E-05 


6.3079E-07 


7.0522E-07 



le; thus, in our case t\ corresponds to t = 13. At ti all 
the exponential terms in equation (3.3) vanish. At ti 
we need to have the terms starting with e~ 2ct vanish; 
this means that ti ~ £i/2, which in our case would 
be 7. Similarly, t% ~ ii/3, which is 4, and so on. 

In Figure 14 (a), (6), (c) and (d), we show the pos- 
terior distributions of A\ , E\ , Ai and c, respectively, 
based on y^ and y-j. The modes of these posterior 
distributions suggest the values of 0.52, 0.029, 0.02 
and 0.4, for E\, A\, Ai and c, respectively. The val- 
ues of A\ and E\ given above are in good agreement 
with the values obtained by a physics based simu- 
lation code. Since the physics based codes are un- 
able to obtain good estimates of A<i and Ei (equiv- 



alently, c), the results on Ai and c obtained by us 
constitute a contribution toward the solution of an 
underlying scientific problem. 

Based on this exercise, plus others that are given in 
Landon (2007), our conclusion therefore is that the 
proposed approach is successfully validated against 
both simulated data as well as the physics code gen- 
erated data. The exercises in Landon (2007) pertain 
to the quark masses of 4 photons and 5 pions. 

6. EXTENDING THE APPROACH 

The approach outlined in Sections 3 and 4 has 
some limitations. The purpose of this section is to 
prescribe strategies for overcoming these. By far, 
the most noteworthy limitation is that the model 
of equation (4.1) restricts attention to a considera- 
tion of the parameters A\, E±, Ai and c, whereas 
the Lattice QCD equations have an infinite number 
of AiS and E^s. The second concern pertains to the 
fact that in Section 5, data associated with the t's 
intermediate to ti and ti are not used in the MCMC 
algorithm. The proposed approach therefore does 
not exploit all the available data yt. Finally, there 
is a question of assuming a constant spacing c of 
the Ei's. What is the effect of unequally spaced E^s 
on inference? Recall that the role played by c is im- 
portant. First, it imparts parsimony by eliminating 
all the EiS save for E\. Second, it gives birth to the 
telescopic series which is central to our approach. It 
turns out that the effect of c is transitionary (it is 
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o 




Fig. 14. Posterior distributions of Ai, E\, A2 and c based on 7/13 and y-j . ( a,) Final posterior of E\. (h) Final posterior 
of A\. (c) Final posterior 0//I2. (d) Final posterior of c. 



a nuisance parameter) and that inferences about E2, 
A3, £3, A4, ... , are possible if we exploit a result ob- 
served in Section 5. 

6.1 Inferences for E2,As and Beyond 

Our ability to extend the approach of Sections 3 
and 4 to the case of E2, A3, E3, A4, . . . , is driven 
by the feature noticed in Section 5.1, that y% does 
not contribute much toward the inferences for A\ 
and Ei, beyond that provided by 7/12. Thus, the ef- 
fect of 1/4, which corresponds to £3 of Table 1, will be 
less so, making it possible for us to do the following: 

Rewrite equation (3.1) as 

00 

(6.1) G {t\.)-A l e- E - t = Y J A n e- Ent , 

n=2 

and let G*(t\-) = G(t\-) - Aie~ Elt , where M and E ± 
are the modes (means) of the posterior distributions 



of Ai and E\ obtained via the likes of Figures 7(b) 
and 8(b). Thus, 

00 

G*(t|.)«£4ie- E » t , 

n=2 

and setting A* n = A n+ \ and E* = E n+ \, for n = 2, 
3, . . . , we have 

00 

(6.2) G*(t\-)*J2 A *n e ~ Kt - 

n=l 

The right-hand side of equation (6.2) parallels the 
right-hand side of equation (3.1), save for the fact 
that A* n and E* replace A n and E n . The mate- 
rial of Sections 3 and 4 now applies, but with the 
caveat that since equation (6.2) is an approxima- 
tion, whereas equation (3.1) is exact, the variance 
of the error terms associated with the former should 
be larger than those associated with the latter. 
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The posterior distributions of £J and A% will be 
the posterior distributions of £2 and A3. The role 
of c as a nuisance parameter is now apparent. The 
posterior distribution of A\ will serve as a revised 
posterior distribution of A2. Indeed, for the MCMC 
runs associated with the treatment of equation (6.2), 
we may sample from the posterior distribution of A2 
to generate the posterior distribution of A*. 

We may continue in the above vein to estimate £3 
and A 4 by defining G**(t\-) = G* {t\-) — A\e~^* 1 , whe- 
re A| and E* are the modes of the posterior distri- 
butions of A* and E±, respectively, and similarly 
with (£4,^15), (E^,Aq), and so on. 

6.2 Using Additional Y t 's 

For enhanced inferences about the parameters A± 
and Ei we may want to use all values of Yt inter- 
mediate to those associated with the labels t\ and ti 
of Tables 1 and 2 and, similarly, with the It's inter- 
mediate to the ones associated with the labels £2 
and £3, and so on. What makes this possible is the 
fact that t\ is the largest value of t for which (A2, E2), 
(A3, £3), . . . , gets annihilated, whereas ti is the lar- 
gest value of t at which (A3, £3), (A4, £4), . . . , gets 
annihilated, and so on. Thus, values of t intermedi- 
ate to t\ and t2 will continue to annihilate (A2,£2), 
(A3, £3), . . . , and those intermediate to ti and t% will 
annihilate (.A3, £3), (A4, £4), . . . , and so on. 

Let uu, y\2, ■ . . , 2/15 denote the Y^s intermediate 
to those associated with the labels t± and £2- Then, 
to incorporate the effect of yn, . . . , j/15 for enhanced 
inference about A\ and E\ , the iterative scheme de- 
scribed in Phase I of Section 4.2 will have to be 
cycled five more times, each cycle involving a use of 
the yij, j = 1, . . . , 5, before proceeding to Phase II, 
wherein the effect of y$ (of Table 1) and 2/7 (of Ta- 
ble 2) comes into play and, similarly, with y2i, the 
single value intermediate to that associated with the 
labels t2 and £3. 

7. SUMMARY AND CONCLUSIONS 

In this paper we have proposed and developed 
a statistical approach for addressing a much discus- 
sed problem in particle physics. Indeed, a problem 
that has spawned several Nobel prizes in Physics. 
The essence of the problem boils down to estimat- 
ing a large (conceptually infinite) number of un- 
known parameters based on a finite number of non- 
linear equations. Statisticians refer to such problems 
as large p — small n. Each equation in our problem 



comprises of the sum of several exponential func- 
tions. 

Previous approaches for addressing this problem 
have been physics based — such as perturbation me- 
thods — and statistics based — such as chi-squared 
goodness of fit, and Empirical Bayes. Physicists have 
found such approaches unsatisfactory, and have cal- 
led for a use of proper Bayesian approaches, thus 
this paper. 

The Bayesian approach proposed by us has been 
facilitated by the fact that by introducing a latent 
parameter, the architecture of the nonlinear equa- 
tions reveals an attractive pattern. This pattern boils 
down to our consideration of a truncated telescopic 
series of equations, each equation being the sum 
of a finite number of exponential functions. Simi- 
lar sets of equations also arise in other arenas of 
science, as mentioned in Section 3.1. The nonlinear 
nature of the equations mandates that our proposed 
approach — which entails stylized proper priors — be 
implemented by a particle style Markov chain Monte 
Carlo (MCMC) approach. Such a procedure turns 
out to be computationally very intensive — about one 
million iterations for making inference about three 
parameters. 

The proposed procedure, when invoked on simu- 
lated data, is able to reproduce the input param- 
eters. This is one way to claim the validity of our 
approach. The procedure, when invoked on some 
real data pertaining to the quark masses of protons 
and pions, is also able to produce results that are 
in agreement with the results produced using alter- 
nate physics based methods. However, the physics 
based methods are able to obtain only partial re- 
sults. By contrast, our approach can produce es- 
timates of as many parameters as is desired — but 
there is no way to validate these against alternate 
approaches or actual numbers, because these are not 
available. 

Future work in this arena will entail enhancements 
to gain computational efficiencies and the choice of 
proper priors that are motivated by a consideration 
of the underlying physics. This means that an un- 
dertaking such as this will call for insights and skills 
that go beyond mathematics, statistics and comput- 
ing. Some appreciation of the underlying physics is 
necessary for, among other things, interest and in- 
spiration! A referee of this paper has made the inter- 
esting suggestion of considering "reference priors." 
This we have been unable to do because, for the 
parameters in question, such priors are not readily 
available. 
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APPENDIX 

Schemata of the 3 Phase MCMC algorithm, which 
can be downloaded at the following: http://www.gwu. 
edu/~stat/irra/Lattice_QCD.htm. 

Phase I 



Specify 
E^andy, 



Generate from 
n(A\E[ a \y,) 



i 

Generate from 



Repeat process to obtain 
(4 2> , ),(4 (3) , ),..., (4 k) , E\ k) ) 



I 



Repeat this whole process m times to 
obtain: 

{a^\e^){a^,e^\..,{a\^e^) 



Go to Phase II 



Fixcatc(°>,4at4 ( f ) 
E, atElWY, at y x and Y 2 at y 2 



Generate A^irom 



T 



Generate A[ from 

n(A\A?,E\\\c^;y x , yi ) 



I 

Generate E[ + ' from 
;r(£j4'UrV»>;^v 2 ) 



I 

Generatec' 1 ' from 

n{c\4\^\Er\y v y 2 ) 
1 



Repeat process to obtain 

(4*\4< 2 *\£< 2 *>, C W ), 
and let this be 

( A(k) A(2k) F (2k) (k)\ 
(■"21 'Al ' a \\ > c l > 



Go to Phase III 



Phase III 



Repeat Phases I and II 
m times using new starting 
values ofii/ ' and c (0) 



We will now have a sample of size 
m from the posterior distribution of 

(A u E v A 2 ,c) 
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